{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# defining the number field K = Q[\\sqrt{5}] with fundamental unit x = (1+\\sqrt{5})/2\n",
    "#\n",
    "K.<x> = NumberField(x^2-x-1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# defining the diagonal quadratic form with coefficients coeff taking vector v as argument\n",
    "#\n",
    "def q(coeff, v):\n",
    "    return sum([c*v^2 for (c,v) in zip(coeff, v)]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# defining the associated bilinear form taking vectors u, v as arguments\n",
    "#\n",
    "def b(coeff, u, v):\n",
    "    return 1/2*(q(coeff, u+v) - q(coeff, u) - q(coeff, v));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# defining the reflection in the hyperplane with normal u\n",
    "#\n",
    "def ref(coeff, u):\n",
    "    n = len(coeff);\n",
    "    v = matrix.identity(n);\n",
    "    return matrix([v[i] - 2*b(coeff, v[i], u)/q(coeff, u)*u for i in range(n)]).transpose();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# computing the order of a matrix m after reduction modulo ideal J = (p) generated by a positive rational prime p in K\n",
    "#\n",
    "def matrix_reduction_order(m, p, max_iter=1E4):\n",
    "    J = K.ideal(p);\n",
    "    k = 1;\n",
    "    l = m.apply_map(J.small_residue);\n",
    "    while not(l==matrix.identity(m.nrows())) and (k < max_iter):\n",
    "        k = k+1;\n",
    "        l = (l*m).apply_map(J.small_residue);\n",
    "    return k;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# dimension = 5\n",
    "#\n",
    "# basis vectors\n",
    "#\n",
    "v = matrix.identity(K, 6);\n",
    "#\n",
    "# coefficients of the quadratic form\n",
    "#\n",
    "q5 = [1-2*x, 1, 1, 1, 1, 1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# outer normals for the polytope P_5\n",
    "#\n",
    "e1 = -v[1] + v[2]; \n",
    "e2 = -v[2] + v[3];\n",
    "e3 = -v[3] + v[4];\n",
    "e4 = -v[4] + v[5];\n",
    "e5 = -v[5];\n",
    "e6 = x*v[0] + (2 + x)*v[1];\n",
    "e7 = x*(v[0] + v[1] + v[2] + v[3]);\n",
    "e8 = (1 + x)*(v[0] + v[1]) + x*(v[2] + v[3] + v[4] + v[5]);\n",
    "e9 = (6 + 9*x)*v[0] + (4 + 7*x)*(v[1] + v[2] + v[3] + v[4]) + (2 + x)*v[5];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# reflections in the hyperplanes of P_5\n",
    "#\n",
    "s1 = ref(q5, e1);\n",
    "s2 = ref(q5, e2);\n",
    "s3 = ref(q5, e3);\n",
    "s4 = ref(q5, e4);\n",
    "s5 = ref(q5, e5);\n",
    "s6 = ref(q5, e6);\n",
    "s7 = ref(q5, e7);\n",
    "s8 = ref(q5, e8);\n",
    "s9 = ref(q5, e9);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# uncomment the following lines to see the reflection matrices\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "#show(\"s1 = \", s1);\n",
    "#show(\"s2 = \", s2);\n",
    "#show(\"s3 = \", s3);\n",
    "#show(\"s4 = \", s4);\n",
    "#show(\"s5 = \", s5);\n",
    "#show(\"s6 = \", s6);\n",
    "#show(\"s7 = \", s7);\n",
    "#show(\"s8 = \", s8);\n",
    "#show(\"s9 = \", s9);\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|delta5|\\phantom{\\verb!x!}\\verb|=| \\left(\\begin{array}{rrrrrr}\n",
       "3 x + 2 & -x - 1 & -x - 1 & -x - 1 & 0 & 0 \\\\\n",
       "0 & 0 & 0 & 0 & 0 & 1 \\\\\n",
       "3 x + 1 & -x & -x - 1 & -x - 1 & 0 & 0 \\\\\n",
       "3 x + 1 & -x - 1 & -x & -x - 1 & 0 & 0 \\\\\n",
       "3 x + 1 & -x - 1 & -x - 1 & -x & 0 & 0 \\\\\n",
       "0 & 0 & 0 & 0 & 1 & 0\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "'delta5 = ' [3*x + 2  -x - 1  -x - 1  -x - 1       0       0]\n",
       "[      0       0       0       0       0       1]\n",
       "[3*x + 1      -x  -x - 1  -x - 1       0       0]\n",
       "[3*x + 1  -x - 1      -x  -x - 1       0       0]\n",
       "[3*x + 1  -x - 1  -x - 1      -x       0       0]\n",
       "[      0       0       0       0       1       0]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "# an infinite order orientation-reversing element in \\Gamma \\cap \\Delta\n",
    "#\n",
    "delta5 = s1*s2*s3*s4*s7;\n",
    "show(\"delta5 = \", delta5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "800"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# the order of the reduction of \\delta_5 modulo (7)\n",
    "#\n",
    "matrix_reduction_order(delta5, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2^5 * 5^2"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# its prime factors\n",
    "#\n",
    "factor(800)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8052"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# the order of the reduction of \\delta_5 modulo (11)\n",
    "#\n",
    "matrix_reduction_order(delta5, 11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2^2 * 3 * 11 * 61"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# its prime factors\n",
    "#\n",
    "factor(8052)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# dimension = 6\n",
    "#\n",
    "# basis vectors\n",
    "#\n",
    "v = matrix.identity(7);\n",
    "#\n",
    "# coefficients of the quadratic form\n",
    "#\n",
    "q6 = [-2*x, 1, 1, 1, 1, 1, 1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# outer normals for the polytope P_6\n",
    "#\n",
    "e1 = -v[1] + v[2];\n",
    "e2 = -v[2] + v[3];\n",
    "e3 = -v[3] + v[4];\n",
    "e4 = -v[4] + v[5];\n",
    "e5 = -v[5] + v[6];\n",
    "e6 = -v[6];\n",
    "e7 = v[0] + x*(v[1] + v[2]);\n",
    "e8 = x*(v[0] + v[1] + v[2] + v[3] + v[4]);\n",
    "e9 = x*v[0] + 2*x*v[1];\n",
    "e10 = (1 + x)*(v[0] + v[1] + v[2]) + x*(v[3] + v[4] + v[5] + v[6]);\n",
    "e11 = (1 + 2*x)*v[0] + (1 + 3*x)*v[1] + (1 + x)*(v[2] + v[3] + v[4]) + x*(v[5] + v[6]);\n",
    "e12 = (1 + 2*x)*v[0] + (2 + 3*x)*v[1] + x*(v[2] + v[3] + v[4] + v[5] + v[6]);\n",
    "e13 = (2 + 2*x)*v[0] + (1 + 2*x)*(v[1] + v[2] + v[3] + v[4] + v[5]) + v[6];\n",
    "e14 = (2 + 3*x)*v[0] + (2 + 4*x)*v[1] + (2 + 2*x)*v[2] + (1 + 2*x)*(v[3] + v[4] + v[5]) + v[6];\n",
    "e15 = (2 + 3*x)*v[0] + (3 + 4*x)*v[1] + (1 + 2*x)*(v[2] + v[3] + v[4]) + 2*x*v[5];\n",
    "e16 = (2 + 4*x)*v[0] + (3 + 6*x)*v[1] + (1 + 2*x)*(v[2] + v[3] + v[4] + v[5]) + v[6];\n",
    "e17 = (3 + 4*x)*v[0] + (2 + 5*x)*v[1] + (2 + 3*x)*(v[2] + v[3] + v[4] + v[5]) + x*v[6];\n",
    "e18 = (4 + 5*x)*v[0] + (4 + 6*x)*v[1] + (2 + 4*x)*(v[2] + v[3] + v[4] + v[5]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# reflections in the hyperplanes of P_6\n",
    "#\n",
    "s1 = ref(q6, e1);\n",
    "s2 = ref(q6, e2);\n",
    "s3 = ref(q6, e3);\n",
    "s4 = ref(q6, e4);\n",
    "s5 = ref(q6, e5);\n",
    "s6 = ref(q6, e6);\n",
    "s7 = ref(q6, e7);\n",
    "s8 = ref(q6, e8);\n",
    "s9 = ref(q6, e9);\n",
    "s10 = ref(q6, e10);\n",
    "s11 = ref(q6, e11);\n",
    "s12 = ref(q6, e12);\n",
    "s13 = ref(q6, e13);\n",
    "s14 = ref(q6, e14);\n",
    "s15 = ref(q6, e15);\n",
    "s16 = ref(q6, e16);\n",
    "s17 = ref(q6, e17);\n",
    "s18 = ref(q6, e18);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# uncomment the following lines to see the reflection matrices\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "#s1 = ref(q6, e1);\n",
    "#s2 = ref(q6, e2);\n",
    "#s3 = ref(q6, e3);\n",
    "#s4 = ref(q6, e4);\n",
    "#s5 = ref(q6, e5);\n",
    "#s6 = ref(q6, e6);\n",
    "#s7 = ref(q6, e7);\n",
    "#s8 = ref(q6, e8);\n",
    "#s9 = ref(q6, e9);\n",
    "#s10 = ref(q6, e10);\n",
    "#s11 = ref(q6, e11);\n",
    "#s12 = ref(q6, e12);\n",
    "#s13 = ref(q6, e13);\n",
    "#s14 = ref(q6, e14);\n",
    "#s15 = ref(q6, e15);\n",
    "#s16 = ref(q6, e16);\n",
    "#s17 = ref(q6, e17);\n",
    "#s18 = ref(q6, e18);\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|delta6|\\phantom{\\verb!x!}\\verb|=| \\left(\\begin{array}{rrrrrrr}\n",
       "814 x + 503 & -279 x - 172 & -179 x - 110 & -178 x - 110 & -178 x - 110 & -178 x - 110 & -4 x - 2 \\\\\n",
       "942 x + 582 & -323 x - 199 & -207 x - 128 & -206 x - 127 & -206 x - 127 & -206 x - 127 & -4 x - 3 \\\\\n",
       "850 x + 526 & -291 x - 180 & -187 x - 115 & -186 x - 115 & -186 x - 115 & -186 x - 115 & -4 x - 3 \\\\\n",
       "420 x + 260 & -144 x - 89 & -92 x - 57 & -92 x - 56 & -92 x - 57 & -92 x - 57 & -2 x - 1 \\\\\n",
       "420 x + 260 & -144 x - 89 & -92 x - 57 & -92 x - 57 & -92 x - 56 & -92 x - 57 & -2 x - 1 \\\\\n",
       "420 x + 260 & -144 x - 89 & -92 x - 57 & -92 x - 57 & -92 x - 57 & -92 x - 56 & -2 x - 1 \\\\\n",
       "64 x + 40 & -22 x - 13 & -14 x - 9 & -14 x - 9 & -14 x - 9 & -14 x - 9 & 0\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "'delta6 = ' [ 814*x + 503 -279*x - 172 -179*x - 110 -178*x - 110 -178*x - 110 -178*x - 110     -4*x - 2]\n",
       "[ 942*x + 582 -323*x - 199 -207*x - 128 -206*x - 127 -206*x - 127 -206*x - 127     -4*x - 3]\n",
       "[ 850*x + 526 -291*x - 180 -187*x - 115 -186*x - 115 -186*x - 115 -186*x - 115     -4*x - 3]\n",
       "[ 420*x + 260  -144*x - 89   -92*x - 57   -92*x - 56   -92*x - 57   -92*x - 57     -2*x - 1]\n",
       "[ 420*x + 260  -144*x - 89   -92*x - 57   -92*x - 57   -92*x - 56   -92*x - 57     -2*x - 1]\n",
       "[ 420*x + 260  -144*x - 89   -92*x - 57   -92*x - 57   -92*x - 57   -92*x - 56     -2*x - 1]\n",
       "[   64*x + 40   -22*x - 13    -14*x - 9    -14*x - 9    -14*x - 9    -14*x - 9            0]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "# an infinite order orientation-reversing element in \\Gamma \\cap \\Delta\n",
    "#\n",
    "delta6 = s7*s13*s18;\n",
    "show(\"delta6 = \", delta6);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# the order of the reduction of \\delta_6 modulo (7)\n",
    "#\n",
    "matrix_reduction_order(delta6, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2^3"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# its prime factors\n",
    "#\n",
    "factor(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# the order of the reduction of \\delta_6 modulo (11)\n",
    "#\n",
    "matrix_reduction_order(delta6, 11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2^2 * 11"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# its prime factors\n",
    "#\n",
    "factor(44)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.9",
   "language": "sage",
   "name": "sagemath"
  },
  "language": "python",
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
